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Variational procedure is developed that yields lowest frequencies of small-amplitude oscillations 
of classical Hamiltonian systems. Genuine Lanczos recursion is generalized to treat related non- 
Hermitian eigenvalue problems. 
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Normal modes £ and frequencies to of small oscillations 
of a classical system near the equilibrium arc determined 
by the secular equation pi 



(1) 



where M and K are N x N symmetric positive definite 
matrices of mass coefficients and spring constants respec- 
tively. In many applications the number N of degrees of 
freedom is large, while only a few lowest frequencies are 
of interest 0. Equation (1) represents a problem more 
complex than a regular symmetric eigenvalue problem, 
unless M or K is diagonal. 

Equation ([!]) can be transformed into the Hamiltonian 
form by introducing canonical momentum r\ — LuAlt;: 



K£ = LOT), T-q = 
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where T = M . Thus, the frequencies of the normal 
modes are the eigenvalues of a 2N x 2N matrix 



T 
K 
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The spectrum of this matrix consists of pairs ±w, since 
(£i ~v) 18 a l so a solution of (||) that corresponds to — u. 
The lowest frequency w m i n is the lowest positive eigen- 
value of the matrix ^ . 

Although the eigenvalues of the matrix (||) are always 
real, the matrix itself is non-Hermitian, unless K = T. 
Therefore, its diagonalization poses a formidable task. 
The major problem is that no general minimum principle 
exists that yields eigenvalues of arbitrary diagonalizable 
non-Hermitian matrices. This does not allow to formu- 
late a variational procedure similar to the Rayleigh-Ritz 
procedure for Hermitian matrices. If K = T, the ma- 
trix (^) is Hermitian, and its positive eigenvalues coincide 
with those of K and T. 

As known from quantum mechanics, the lowest eigen- 
value e m i n of a Hermitian matrix H can be obtained from 
the minimum principle 
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The minimum is to be searched over all vectors ^. The 
Ritz variational procedure is an approximation when the 



set {ip} in (Q) is restricted to some subspace K, of dimen- 
sion n < N. The best approximation to e m i n in the sense 
of ([|) is obtained as the lowest eigenvalue of the n x n 
Rayleigh matrix H, obtained by projection of H onto IC. 

The special paired structure of the matrix (j^) makes it 
possible to generalize (||) such as to yield w m i„. In fact, 
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The minimum is to be searched over all possible phase 
space configurations {£,,rf\. Before providing the proof 
to this equation let me point out some of its features. 

First, it states that w m i n is the minimum harmonic part 
of the total energy, (£K£)/2 + (r]Tr])/2, over the phase 
space configurations normalized by (£77) = 1. Since K 
and T are both positive definite, the right-hand side is 
strictly positive and so is w m in- Second, equation (||) is 
symmetric in K and T, according to the nature of the 
problem. When K = T the minimum is achieved at 
£ = 7], and (||) becomes the same as (f|). 

Note that the functional in (||) has no maximum, since 
the denominator can be made arbitrarily small. The 
global minimum, however, always exists. This is not ob- 
vious, since a set of all pairs of vectors with = 1 is 
not compact. Indeed, say, any vector orthogonal to 77 can 
be added to £, making |£| arbitrarily large. However, the 
functional in (|5|) grows indefinitely in this case, such that 
the global minimum is achieved at finite |£| and |7y|. 

Variation of (||) with respect to £ and 77 yields equa- 
tions (||). Thus, the solutions of (||) are the stationary 
points of (H). The global minimum (|^), therefore, gives 
indeed w m i n . The singularity in the denominator poses 
no problem, since it corresponds to infinitely large values 
of functional, while near the minimum it is analytic. 

Minimum principle (|J) can, in fact, be obtained from 
the Thouless minimum principle ||, derived for non- 
Hermitian matrices that appear in random phase approx- 
imation (RPA). Equation (||) transforms into the Thou- 
less minimum principle by substitution: A = (K + T)/2, 
B = (K- T)/2, x = (£ + 77)/2, and y = (£ - r?)/2. 

Variational procedure similar to the Rayleigh-Ritz pro- 
cedure can be formulated if coordinates £ and momenta 
77 in (H) are restricted to some subspaces U and V of di- 
mension n, respectively. 
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Let and {?7i} be two sets of vectors that span U 
and V, such that = Expanding £ = 

77 = v i r li an< i varying (@) with respect to Ui and Vi, we 
find the latter to obey a 2n x 2n eigenvalue equation 

(!:)(:)-(:)■ 

with ify = (S,iK^j) and = (ijiTrjj). Equation (||) 
generalizes Hermitian Rayleigh-Ritz eigenvalue equation 
for H. It has 2n solutions ±ui, the lowest positive of 
which gives the best approximation to w m i n in the sense 
of equation (||) . 

Krylov subspace Q for the matrix ^ can be con- 
structed by acting with it many times on an arbitrary 
vector (£1,771): 

($)■ (30- (£?£)■■■• < 7 > 

The subspace that spans first n vectors of this sequence 
has the property of approximating an invariant subspace 
of (||). Thus, it is natural to expand approximation to 
an eigenvector of (||) as a linear combination of these 
vectors. In other words, the natural choice for the sub- 
spaces IA and V for the variational procedure described 
above are the subspaces Li n and V n that span the upper 
and lower components of first n vectors of (0) . 

In order to implement the variational procedure, it 
is necessary to construct a biorthogonal basis {^i,rji}, 
i = 1, ■■■ 1 n in lA n and V n and compute matrix elements of 
K and T . Both tasks can be performed simultaneously 
using the following recursion: 

6+1 = Pr+i (TVi - a£i ~ (8a) 
Vi+l = $i+l( K €i - 7iVi - 8iVi-l)- ( 8b ) 

The four coefficients ctj, fy, 7$, and Si are to be cho- 
sen at each step i such as to make orthogonal to r\i 
and rji-i, and rji+i — orthogonal to £j and This 
appears to be enough to ensure global biorthogonality 
(kiVj) = <%■ 

Indeed, assume biorthogonality to hold up to step i. 
Multiplying (|8a| ) by r\j, j < i — 1, we have (%£i+i) oc 
(rjjTrji) = (rjiTrjj) = due to Hermiticity of T and 
the fact that Tr\j is a linear combination of all with 
k < j + 1 < i. Thus, the biorthogonality also holds for 
the step i + 1. 

Multiplying (|a|) by r?i_i, 77*, and rji + i and using 
biorthogonality, we get (£i?7i) = 1, Ku = on, and 
Ki^-i = Ki-\,i = Pi- Similarly, T u = 7, and Ti,i-i = 
Ti_i.i = Si. All other matrix elements of K and T vanish. 

The recursion (^J) is a straightforward generalization of 
the Hermitian Lanczos recursion [||,Q 

= ^^(Hipi - a % ipi - fa-ipi-i) (9) 



applicable to any Hermitian matrix iJ. When K = T 
and £1 = 771 both equations (||) coincide with each other 
and with equation ([)]), up to the notation. 

As in the case of the Hermitian Lanczos algorithm, 
several lowest frequencies can be found one by one by 
projecting the £— and 77-components of converged eigen- 
vectors out of V n and U n subspaces respectively. 

The method was tested on a set of large sparse ran- 
dom matrices of the form (|J). Symmetric matrices T 
and K were generated to have an average of 40 randomly 
distributed and randomly positioned matrix elements in 
each row. Both K and T were shifted by an appropri- 
ate constant to ensure positive definiteness. Figure 1 
demonstrates the convergence results for a matrix of the 
size 2N = 200000. 
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FIG. 1. Convergence of the generalized Lanczos algorithm 
for a random matrix of the form (^) and size 2N = 200000 

For smaller matrices up to 2A?~ = 2000, where it was 
possible to obtain all eigenvalues with regular methods, 
the present method has converged to the true lowest fre- 
quency in all instances. 

In conclusion, the method is proposed that generalizes 
Rayleigh-Ritz variational procedure and Lanczos recur- 
sion to the case of non-Hermitian matrices of the form 
, that determine normal modes and frequencies of 
small-amplitude oscillations of Hamiltonian systems. 

Equations (^) have numerous applications beyond 
purely mechanical problems. Schrocdinger equation in 
non-orthogonal basis represents a generalized symmetric 
eigenvalue problem similar to ([!]). RPA and other time- 
dependent techniques in nuclear physics and quantum 
chemistry lead to the equations similar to (||) At 
last, eigenvectors of so-called Hamiltonian matrices, to 
which (||) is a special case, solve the nonlinear algebraic 
Riccati equation which appears in the theory of stability 
and optimal control Q. 

I would like to acknowledge numerous enlightening dis- 
cussions with Vladimir Chernyak during my appointment 
at the University of Rochester. 
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